BB103E 


00  IN 


t>  u 


lUL 


TM  No.  861032 


NAVAL  UNDERWATER  SYSTEMS  CENTER 
NEW  LONDON  LABORATORY 
NEW  LONDON,  CONNECTICUT  06320 


Technical  Memorandum 


nr ;  rnriipr 

LflLiiOw 


ONLY 


COMPUTATION  OF  COMPLEX  AIRY  FUNCTIONS 


Date:  19  March  1986 


Prepared  . 

'  ■»»*. _ i  .«*  » .  * 


Marvin cJ.  Goldstein 
Surface  Ship  Sonar  Dept 


L.ERENC  ONLY 


Approved  for  Public  Release:  Distribution  Unlimited 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

19  MAR  1986 


2.  REPORT  TYPE 

Technical  Memo 


3.  DATES  COVERED 

19-03-1986  to  19-03-1986 


4.  TITLE  AND  SUBTITLE 

Computation  of  Complex  Airy  Functions 


6.  AUTHOR(S) 

Marvin  Goldstein 


5a.  CONTRACT  NUMBER 
5b.  GRANT  NUMBER 
5c.  PROGRAM  ELEMENT  NUMBER 
5d.  PROJECT  NUMBER 

W65000 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Underwater  Systems  Center, New  London, CT, 06320 


5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

TM  No.  861032 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES )  10.  SPONSOR/MONITOR' S  ACRONYM(S) 

1 1.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

NUWC2015 

14.  ABSTRACT 

In  this  paper,  we  discuss  an  efficient  way  to  evaluate  scaled  complex  Airy  functions  by  asymptotic  series. 

By  obtaining  an  a-priori  estimate  of  the  number  of  monotonically  decreasing  terms  of  the  asymptotic  series 
for  exp  (zA1.5/1.5)Ai(z)  that  can  be  summed  without  underflowing  the  unit  round-off  of  the  computer 
arithmetic,  we  avoid  unnecessary  computations  during  the  summation.  Using  this  method,  we  develop  an 
efficient  computational  routine  for  obtaining  numerically  linearly  independent  Airy  functions  whose 
Wronskian  is  stable  throughout  the  complex  plane.  The  scaled  complex  Airy  function  of  the  second  kind 
exp  (-zA1.5)Bi(z)  may  be  obtained  with  well-known  connecting  formulas  from  the  numerically  independent 
computations.  However,  the  Wronskian  of  the  Airy  functions  of  the  first  and  second  kind  is  not  stable  in 
the  sector  180??>  abs(  arg  z  )  >  60??,  when  abs(z)  is  sufficiently  large. 

15.  SUBJECT  TERMS 

Airy  Function;  Wrongskian 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

25 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


REFERENCE  ONLY 


TM  No.  861032 


NAVAL  UNDERWATER  SYSTEMS  CENTER 
NEW  LONDON  LABORATORY 
NEW  LONDON,  CONNECTICUT  06320 


Technical  Memorandum 


COMPUTATION  OF  COMPLEX  AIRY  FUNCTIONS 


Date:  19  March  1986 


Prepared  (?.  . 

Marvin  cJ.  'Goldstein 
Surface  Ship  Sonar  Dept 


Approved  for  Public  Release:  Distribution  Unlimited 


TM  861032 


ABSTRACT 


In  this  paper,  we  discuss  an  efficient  way  to  evaluate  scaled  complex 
Airy  functions  by  asymptotic  series.  By  obtaining  an  a-priori  estimate  of 
the  number  of  monotonical ly  decreasing  terms  of  the  asymptotic  series  for 
exp  <z,s/'l ,  5>Ai  (z)  that  can  be  summed  without  underflowing  the 
unit  round-off  of  the  computer  arithmetic,  we  avoid  unnecessary 

computations  during  the  summation.  Using  this  method,  we  develop  an 
efficient  computational  routine  for  obtaining  numerically  linearly 
independent  Airy  functions  whose  Wronskian  is  stable  throughout  the  complex 
plane.  The  scaled  complex  Airy  function  of  the  second  kind 

exp  <—z’-s/'l .  5)Bi  (z)  may  be  obtained  with  well-known  connecting 
formulas  from  the  numerically  independent  computations.  However,  the 
Wronskian  of  the  Airy  functions  of  the  first  and  second  kind  is  not  stable 
in  the  sector  180e >  abs<  arg  z  )  >  60°,  when  abs(z)  is 
sufficiently  large. 


ADMINISTRATIVE  INFORMATION 


This  memorandum  was  prepared  under  Job  Order  No.  W65000,  EVA  Prog. 
Modeling  Efforts,  Principal  Investigator,  H.  Weinberg  (Code  3332).  The 
author  of  this  memorandum  is  located  at  the  New  London  Laboratory,  Naval 
Underwater  Systems  Center,  New  London,  Ct.  06320. 
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INTRODUCTION 

Typically  the  asymptotic  expansion  of  the  scaled  Airy  function 

exp  (zta)Ai  (z)  <zta=z,s/'l ,  5)  is  truncated  to  obtain  a  function 
that  approximates  exp(zta) Ai  (z)  with  arbitrary  accuracy  for 
sufficiently  large  values  of  abs<z>.  On  the  other  hand,  for  any 
specific  value  of  z  there  is  a  limit  to  the  accuracy  with  which  we 
can  compute  exp  (zta) Ai  (z)  by  truncating  its  asymptotic  expansion,  and 
the  accuracy  cannot  always  be  improved  by  taking  additional  terms  of  the 
asymptotic  expansion.  Accuracy  may  be  improved  by  adding  additional  terms 
as  long  as  the  magnitude  of  the  terms  continue  to  decrease  without 
underflowing  the  computer’s  unit  round-off  error.  Although  Thomson  [1]  has 
implemented  a  complex  Airy  function  computer  code,  the  code  runs  slowly 
because  comparison  tests  are  performed  to  assure  monotonicity  and  avoid  the 
accumulation  of  terms  that  underflow  the  smallest  positive  floating-point 
computer  number  instead  of  the  computer’s  larger  unit  round-off  error. 

In  this  paper,  we  discuss  an  efficient  way  to  evaluate  scaled  complex 
Airy  functions  by  asymptotic  series.  By  obtaining  an  a-priori  estimate  of 
the  number  of  monotonically  decreasing  terms  of  the  asymptotic  series  for 
exp  (zta) Ai  (z)  that  can  be  summed  without  underflowing  the  unit 
round-off  of  the  computer  arithmetic,  we  avoid  unnecessary  computations 
during  the  summation.  Using  this  method,  we  develop  an  efficient 
computational  routine  for  obtaining  numerically  linearly  independent  Airy 
functions  whose  Wronskian  is  stable  throughout  the  complex  plane.  The 
scaled  complex  Airy  function  of  the  second  kind  exp(-zta)Bi  <z)  may  be 
obtained  with  well-known  connecting  formulas  from  the  numerically 
independent  computations.  However,  the  Vronskian  of  the  Airy  functions  of 
the  first  and  second  kind  is  not  stable  in  the  sector  180e >  abs(  arg  z  ) 

>  60°,  when  abs(z)  is  sufficiently  large. 
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ASYMPTOTIC  ANALYSIS 


Consider  approximating  the  scaled  Airy  function  exp  (zta)Ai  (z) 
for  values  of  z  of  large  magnitude  by  evaluating  a  partial  sum  of  the 
asymptotic  series  C2I 

exp(zta)Ai  (z>  ~  (n-°az-°3S/2)  2  (-l)kckzta~k  (1) 

A-o 

where 


zta  =  z,s/1.5,  abs<  arg  z  )  <  rr 
Co  =  1,  ck+,/ck  =  k/2  +  5<k+l>-’/'72 .  <2> 

Writing 

Gk  =  (-1 )  kckzta  '* 

we  see  that  -fche  terms  of  the  sum  in  <1>  can  be  generated  recursively: 

<?*.,  =  -  (ck+,/'ck)zta~'Gk  (k  =  0,1,2,  ...)  (3) 

and,  a  fortiori, 

7^  < Cj ♦/ 


abs(Gk+,)  =  abs(zta~k~’) 


=  2~k}zl abs  (zta~k~')5/72  =  abs(Gk+i).  (4) 

Let 

✓7  -  / 

S„.,(z)  =  ][_Gk  .  (5) 

By  definition  of  an  asymptotic  series  [31,  as  z  ->  °o 
lim  ztan(2n  sz3sexp  (zta)  Ai  <z)~S„.,  (z)  )  =  (~l)nc„  . 

Therefore,  for  values  of  z  of  sufficiently  large  magnitude 
2noaz03Sexp  (zta) Ai  (z)  -  S„-t(z)  =  G„  .  (6) 

For  a  given  value  of  z,  if  G„  is  the  first  term  in  series  S^(z) 
for  which 


abs(G„ .,)  >  abs  (G„)  (7) 

then  the  magnitude  of  the  error  in  S„-,  (z)  is  ’’minimized” ,  since 
G„  will  have  the  smallest  magnitude  of  all  the  terms  in  the 
series.  By  <3>,  condition  <7>  is  equivalent  to 
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c„„/c„  >  m  =  abs  (zta) .  (8) 

Therefore,  the  first  value  of  n  for  which  <7>  holds  is 
n  =  [  (9m3  +  9m  +  l)°*/3  +  m  -  0.5J~[2m  -  (14.  4m)-'J  (9) 

where  Cxi  is  the  smallest  integer  greater  than  or  equal  to  x. 
However,  if  m  is  sufficiently  large,  the  computation  of  G„  will 
cause  computer  underflow  well  before  k  assumes  the  value  of  n  given 
in  <9>,  and  the  evaluation  of  <5>  out  to  then  becomes 

impractical.  In  this  case,  <5>  could  be  summed  until  abs(G„) 
underflows  the  smallest  positive  computer  number  s.  In  14],  it  was 
shown  how  to  accomplish  this  without  executing  comparison  tests  on 

G„,  resulting  in  a  40%  improvement  in  execution  efficiency. 

However,  an  additional  50%  improvement  in  execution  efficiency  is  passible 
by  considering  the  asymptotic  behavior  of  the  relative  error  in  the 
approximation  Sn.,  (z)  to 

2nosz°-3Sexp(zta)Ai(z).  <10> 

Since 

S„.,(z)  -  2ir°  sz°  3Sexp  (zta)Ai  (z)  ->  0  as  z  ->  (11) 

it  follows  that 

2/?°  *ze  ssexp(zta)Ai  (z)  ->  1  as  z  ->  **<=  (12) 

and,  a  fortiori,  for  the  relative  error  in  S„-t(z)  we  have 
S„.,  (z)  -  2n0  Sz0  2Sexp(zta)Ai(z) 

_ : -  - >  (S„.,(z)  -  1)/1 

2tt°  *zo  ssexp  (zta) Ai  (z) 

- >  0/1  (13) 

as  z  ->*o.  But  the  difference  between  a  function  and  any  partial  sum  of 
its  asymptotic  series  is  of  the  order  of  the  first  neglected  term  of  the 
series  as  z  -><*> .  Therefore,  <13>  suggests  that  for  sufficiently  large 
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values  of  ahs(z),  we  should  be  able  to  approximate  <10>  almost  to  the 

relative  precision  of  the  computer  arithmetic  by  summing  <5>  until 

abs(Gk)  is  smaller  than  the  computer’s  unit  round-off  error  M 

instead  of  the  computer’s  smallest  positive  number  s.  For  a 

practical  computer  evaluation  of  S„.,(z)  that  avoids  comparison 

tests  on  <?*,  we  develop  an  a-priori  estimate  that  depends  on  z 

and  M  of  the  number  of  monotonically  decreasing  terms  of  S„.,(z) 

that  can  be  summed  without  underflowing  the  computer’s  unit  round-off 

errors  ,  To  do  this,  we  proceed  in  essentially  the  same  way  we  did  in  [4], 

replacing  the  smallest  positive  computer  number  s  with  the  computer 

unit  round-off  error//, 

An  estimate  of  the  largest  value  of  m=abs  (zta)  for  which  the 
magnitude  of  the  terms  of  the  partial  sum  S„(z)  is  monotonically 
decreasing  and 

n 

M<  abs  <G„)  <  abs(Gn *,  ;  =  #[7(1  +  (j+l)-'/(7.  2j  )  >  (14) 

can  be  obtained  by  setting  approximation  <4>  for  abs  (G„.,)  equal 

to  /<  with  n=2m,  applying  Stirling’s  asymptotic  formula  for  factorials 

and  solving  the  resulting  equation  for  m.  This  gives  m*  as  the 

largest  value  of  m,  where  zo*  is  the  limit  of  the  rapidly 

convergent  sequence 

+  me,  mo=-ln(7.  2a*h~c  *)/2  (15) 

In  other  words,  for  any  value  of  z  for  which  m<m*, 
the  theoretical  ’’minimal”  error  in  Sn~,  (z)  (namely  abs(G„) 
with  n=2m)  is  computationally  attainable,  because  it  exceeds*/.  On  the 
other  hand,  for  any  value  of  z  for  which  m >m*,  the  theoretical 
’’minimal”  error  in  Sn.,(z)  will  underflow*/  or  be  close  to*/, 
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since  it  satisfies 

•5 

abs(G„(m))  <  abs(Gk(m))  <k  =  2m*) 

.  a 

=  abs  <Gk  <m)  ) 

=  abs(Gk*,< m*))(m/m*)~k  ,  by  (4) 

=  M  (m/m*)  '*  ->  C  as  m/m*  ->o£>.  <  16) 

Therefore,  for  values  of  z  for  which  m>m*,  we  can  approximate 

<10>  almost  to  the  relative  precision  of  the  computer  arithmetic  by  summing 

no  more  than  k=2m*  terms  of  #S,-»  <z) j  in  fact,  we 

evaluate  SCkj-p-i  (z)  for  a  positive  integer  value  of 

p<[kJ,  in  order  to  avoid  computing  terms  of  SCkJ.,  (z)  that 

underflow  m.  when  m/m*  is  large.  ' 

To  select  StkJ-p (z)  when  m>m*,  we  determine 
the  index  [ kJ -p  <=[ k-pJ )  of  the  error  term 

M-p-* 

abs  (GIkJ.p  <m) )  =  a  bs  <GCk].p  (m) )  II  (  1  +  (J  +  l)-'/(7.2j)  ) 

for  which 

A  A 

/< <  abs(G[k].<p*t)  <m)  )  >  abs  (G[k]-P(m)  )  < u/  (17) 

(p  <  EkJ,  k  =  2m*) 

in  the  following  manner. 


i 

then  eq.(15)  gives  a*=19.35  (abs(z)=9.44); 

therefore,  for  absfz)  >  9.44,  no  sore  than  39  teras  of  S„-i(z)  are 
required  to  approbate  (!0)  to  the  relative  precision  of  the  coaputer 
arithaetic.  The  nuaber  of  required  teras  falls  off  at  a  soderate  rate  as 
a/a*  aoves  away  froa  !;  e.g.  for  a/a*=2,  5.7  (abs(z)=15 ,  30), 
the  nuaber  of  teras  required  is  15,  10,  respectively.  A  list  of  the 
required  nuaber  of  teras  [2a*-p]  for  corresponding  values  of  a  and 
a/s*  appears  in  Table  1 . 
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Obtain  a  table  of  values  of  m,  namely  m?,  for  positive 
integer  values  of  p<[kJ,  which  satisfy  the  equation 

A 

abs  (Gtti-P(m)  )  =M.  (18) 

For  a  given  value  of  m>m*,  say  m  =  m,  If  mp<m<mp*,, 
then 

*  =  abs(G,k].p(mp))  >  a bs  (Gck}.p  (m) )  (19) 

and 

A  A  ^ 

M  =  abs  (GCi,]-(p,,,  (mp*, ))  <  abs  (Giu-r-t  (m) )  (20) 

since 

A 

abs  (G[kJ.p(m)  )  >/U,  if  m  <  nip 

<JU  ,  if  m  >  m,  . 

Therefore,  if  n^<S<mp*i  we  evaluate 

A  _ 

St nj-(p-n)  (z) ,  since  Gcm-p  (m)  satisfies  <17> ,  by 
<19>  and  <20>  .  For  example,  assuming  AC  =2’*°,  if  a  given  value 

of  m>m*,  say  m=19.6,  lies  between  the  values  of  m  for  p=5,6  in 
Table  1,  then  take  [2m*-53=34  for  C k-pl . 
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p 

m 

m/m* 

[2m* 

0 

19.35258 

1 

39 

1 

19.36194 

1. 000484 

38 

2 

19.38579 

1. 001716 

37 

3 

19.42577 

1. 003782 

36 

4 

19.48381 

1 . 006781 

35 

5 

19.5621 

1. 010827 

34 

6 

19.66328 

1. 016055 

*  33 

7 

19.79037 

1. 022622 

32 

8 

19.94698 

1. 030714 

31 

9 

20. 1374 

1. 040554 

30 

10 

20.36676 

1.052406 

29 

11 

20.64125 

1.066589 

28 

12 

20.96842 

1. 083495 

27 

13 

21.35754 

1. 103602 

26 

14 

21.8201 

1. 127504 

25 

15 

22.37055 

1. 155947 

24 

16 

23.02719 

1. 189877 

23 

17 

23.81359 

1.230513 

22 

18 

24.76054 

1.279444 

21 

19 

25.90883 

1.338779 

20 

20 

27.31371 

1.411373 

19 

21 

29.05135 

1.501162 

18 

22 

31.22937 

1.613706 

17 

23 

34.00401 

1.757079 

16 

24 

37.60926 

1 . 943372 

15 

25 

42.4084 

2. 191357 

14 

26 

48.99014 

2.531453 

13 

27 

58.35861 

3. 015547 

12 

28 

72.336 

3.737796 

11 

29 

94.49151 

4.882632 

10 

30 

132.5262 

6.847986 

9 

31 

205.2731 

10. 60702 

8 

32 

367.2325 

18.97589 

7 

33 

818.2973 

42.28363 

6 

34 

2605.586 

134.6377 

5 

35 

15653.95 

808. 882 

4 

36 

342085.9 

17676.5 

3 

37 

2. 000795E+08 

1. 033865E+07 

2 

38 

8. 006342E+16 

4. 137093E+15 

1 

Number  of  terms  needed  to  obtain  exp(zta)  Ai(z) 
to  a  relative  precision  of^=  2"^  for  values  of 
m>m*  =  19.35  is  [2m*  -  p] 

TABLE  1 
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Analyzing  the  asymptotic  series  for  exp  (zta)Ai '  (z)  : 

exp(zta)Ai'  (z) - (n-°*z0Xi/2)  X  (~l)kdkzta-k  (21) 

k=o 

where 

zta  =  z's/1.5,  abs  (  arg  z  )  <  n 
d0  =  1,  dk.,/dk  =  k/2  -  7<k+l)-’/72  (22) 

in  the  same  way  we  analyzed  <1>  produces  similar  results.  Writing 

Gk  =  (-1 )  kd„zta  ** 

we  see  that  the  terms  of  the  sum  in  <21 )  can  be  generated  recursively: 
Gkt,  =  ~(dk.,/dk)zta-’Gk  (k  =  0,1,2,  ...)  (23) 


and,  a  fortiori, 

k 

abs(Gk+,)  =  abs  (zta~k't)  (7/72)  JT  (dik,/dt) 

=  2~kk!  abs  (zta~k~’)  7/72  =  abs  (Gk.,)  (24) 

Let 

/>-/ 

S„.,(z)  =  X  Gk  (25) 

k-0 

For  a  given  value  of  z,  the  magnitude  of  the  theoretical  error  in 
S„.,(z)  (namely  abs(G„))  is  minimized  when 


n  -  [  2m  +  7/(12m)J 


An  estimate  of  the  largest  value  of  m  for  which  the  minimal  error  satisfies: 
n-f 

abs  (G„)  =  z/,TT(  1  ~  14  (J  +  1 )  ~'/(72j )  )  (**‘  =  7<U/5) 

J- 1 

can  be  obtained  by  setting  approximation (24)  for  abs(G„)  equal 
to^'with  n-l=2m,  applying  Stirling’s  asymptotic  formula  for 

factorials  and  solving  the  resulting  equation  for  m.  This  gives 

■®*  as  the  largest  value  of  m,  where  m*  is  the  limit 
of  the  rapidly  convergent  sequence 
mJ=-ln(mJ-,)/4  +  mo,  mo=-ln  (7 .  2jurr0  S) /2 


12 


TM  861032 


Therefore ,  for  values  of  x  for  which  m>m*,  we  can 
approximate 

-  <2*1*  *x-°  *s)exp<zta)Ai  '  (z) .  <27 > 

to  the  relative  precisian  of  the  computer  arithmetic  by  summing  no  more 
than  k=2m*+l  terms  of  S„-i<z>;  in  fact,  for  a  positive 
integer  value  of  p<[kJ,  we  evaluate  SCkJ.p.,  (x) ,  where 
[  k] -p  is  the  index  of  the  error  term 

a  IkhP't 

ahs  <G„3.P  (m>>  =  abs  <GCkJ.„  <m) )  ff  <  1  -  14  (J  +  l)-’/<72j  )  ) 

J*' 

that  satisfies 

/  +  *  / 

^  abs  ( j-c p+i j  (m)  )  )  abs  (Gckj-P  (m)  )  C  (28) 

(p  <  i  k  J ,  k  =  2xo*  +  1) , 

The  index  CkJ—p  is  determined  in  the  following  manner. 

Since  <4t^=7.&'/5,  the  table  of  values  of  m  for  which  approximation <24 > 
reduces  to  A*is  identical  to  Table  1  far  <et  =2~*°.  Therefore,  for 
a  given  value  of  m>m*,  we  determine  CkJ-p 

<=[  2m*3  -p+1)  by  locating  where  m  falls  in  Table  1.  For  example,  if 
a  given  value  of  m>m *,  say  xa=19.  6,  lies  between  p=5,6  in  Table 
1 ,  take  [ 2m*— 4-3  =35  for  [kJ-p  ( =[  2m* 3  -p+1  >  . 

COMPUTATIONAL  STABILITY  OF  THE  VRONSKI AN 

As  pointed  out  in  1 5]  ,  in  the  sector  60°<  abs(arg  x)  <180° 
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Al(z)~  +iBi(z),  Ai'(z)~  +  iBi(z)  (abs(z)  large)  (29) 

Therefore,  although  theoretically  the  Wronskian  I V(Al(z),  Bl(z)) 
satisf les 

V(Al(z),  Bi(z))  =  n",  (30) 

it  is  numerically  unstable  In  this  sector  when  abs(z)  is  large 
because  of  severe  cancellation  of  significant  digits;  in  other  words, 
computationally  V(A1  (z) ,  B1  (z)  )  may  be  close  to  zero.  For  example,  if 
abs  (z)=15  and  arg  z=150°,  then  (to  at  least  20  significant 
digits) 

Beal  Part  Imaginary  Part 

Al(z):  -3010621074.950100596341  112308057590.  0163025659  . 

B1  (z) :  -112308057590.  0163025659  -301 0621 074.  9501 00596341 

Ai '  fz;  :  422189739015.  1793969374  -99699231497.  52144586569 

B1‘  (z):  99699231497.52144586569  422189739015.1793969374 

so  that 

Ai(z)  =  -iBl(z),  Ai'  (z)  =  -iBi'  (z)  (31) 

and  V(Ai(z),  Bl(z))=0  on  computer  aithmetic  hardware  having  20 
significant  decimal  digits  or  less.  On  the  other  hand,  this  instability 
does  not  occur  for  the  following  numerically  linearly  independent  Airy 
functions 

Ai(z),  Al(z  exp  (±2ni/3)  )  (32) 

where 

V(Ai(z),  Ai  (z  exp(2ni/3))  =  0.  5n~'exp(-nl/6)  (Im(z)  <  0) 

W(Ai(z),  Ai  (z  exp(-2ni/3))  =  0.  5n~,exp(vi/6)  (Im(z)  >  0) 

14 
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COMPUTER  PROGRAM 

For  m  >  19.35  (abs(z)  >  9.44),  the  complex  Airy  function 

subprogram  computes  the  scaled  numerically  linearly  independent  Airy 
functions: 

exp(zta)  Ai(z),  exp(-zta)  Ai  (z  exp(+2nl/3))  (33) 

and  their  scaled  derivatives 

exp(zta)  Ai’(z),  exp(-zta)  Al  ’  (z  exp  (+2ni/3 )  )exp  (p2ni/3)  (34) 

by  evaluating  partial  sums  of  their  asymptotic  expansions  <1,  21)  to  the 
working  precision  of  a  computer  with  unit  round-off  error 
M=2-*°,  avoiding  the  accumulation  of  terms  that  underflows.  In 
evaluating 

exp(-zta)  Ai  (z  exp  (~2ni/3)  ) ,  exp(-zta)  Ai  '  (z  exp  (+2ni/3)  )  exp  (-2/T1/3) 
we  rotate  z  120°  counter — clockwise  if  its  imaginary  part  is 
negative;  otherwise  z  is  rotated  120°  clockwise  away  from 
180°. 

For  7.45  <  m  <  19.35  (5  <  abs(z)  <  9.44),  we  evaluate  only 

[2m+ll  terms  of  the  asymptotic  expansions  <1,  21),  so  that  the 

largest  error  occurs  for  m=7.  45 

For  0  <  abs(z)  <  5,  we  employ  the  ascending  power  series  given 
in  C21  to  compute  the  Airy  functions  of  the  first  and  second  kind  from 
which  we  obtain 

exp(0)  Ai  (z  exp  (+2ni/3) )  =  0.5exp(+ni/5)  (Bl(z)  +  iAi  (z) ) 

exp(O)  Ai  '  (z  exp  (+2nl/3)  )exp(+2nl/3)  =  0.  5exp(-nl/'6)  (Bi'(z)  +  lAi’(z)) 
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The  scaled  complex  Airy  function  of  the  second  kind: 
exp(-zta)Bi  (z) ,  exp(-zta)'Bi ’  (z)  (abs(  arg  z  )  <  60°)  (35) 

exp  (zta)Bl  (z) ,  exp(zta)Bl'  (z)  (60°  <  abs(  arg  z  )  <  180 °) 

may  be  obtained  in  the  asymptotic  region  (abs  (z)  >5)  from  the 
numerically  linearly  independent  scaled  computations  by  utilizing  the 
connection  formulas: 

Bi(z)  =  2exp(+rd/6)  Ai  (z  exp(+2ni/3))  p  lAl(z)  (36) 

Bi  '  (z)  =  2exp  (+ni/6)  Ai’  (z  exp(+2rti/3)  )exp(+2ni/3)  +  1A1  ’  (z) 

If  abs(  arg  z  )<60e,  then  Real  (zta)  >0,  so  we  compute 

exp  (-zta)  Bi  (z )  =  2exp(+rti/6)  (  exp  (-zta)  Ai  (z  exp  (+2ni/'3))  ) 

^ iexp(-2zta )  (  exp  (zta)  Ai  (z)  )  (37) 

exp  (-zta)Bi  ’  (z)  =  2exp(+5ni/6)  (  exp  (-zta)  Ai’  (z  exp  (+2ni/3)  )  ) 

^  iexp(-2zta)  (  exp(zta) Ai’  (z)  ), 

where  the  upper  sign  is  selected  if  lm(z)<0,  and  the  lower  sign,  if 
lm(z)>0 .  If  60°<abs(  arg  z  )  <180°,  then  Real  (zta)  <0, 
so  we  compute 

exp  (zta)  Bi  (z)  =  2exp(±ni/6)  exp(2zta)  (  exp(-zta)Ai(z  exp(+2ni/3))  ) 

+.  i  <  exp  (zta) Ai  (z)  )  (38) 

exp  (zta)  Bi’  (z)  =  2exp(+5ni/6)  exp(2zta)  (  exp  (-zta)  Ai  ’  (z  exp(+2ni/'3 ))  ) 

+  i(  exp  (zta)  Ai  ’  (z)  ), 

where  again  the  upper  sign  is  selected  if  lm(z)<0,  and  the  lower 
sign,  if  lm(z)>0. 
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***************************************************************** 

IMPROVED  COMPLEX  DOUBLE  PRECISION  AIRY  FUNCT ION . . . M .  GOLDSTEIN 


GIVEN  Z*  IF  INDEP  >  0.  CLAIRY  RETURNS  THE  NUMERICALLY  LINEARLY 
INDEPENDENT  AIRY  FUNCTIONS: 

B I  =  EXP(-ZTAM)  *  A 1 C  Z*EXP(-/+ \ 2PI /3)  ) 

BIP  =  EXP(-ZTAM)  *  A1(  Z*EXP ( -/+ 1 2PI /3 )  )  •  EXP( -/+ 1 2PI /3 ) 

(  AI,  A I P  )  =  EXP(ZTAM)  •  (  AI(Z),  A1(Z)  } 

(  ZTAM  =  2  *  SORT C  Z )  *  *3  /  3  ) 

FOR  THE  COMPUTATION  OF  BI  AND  BIP,  Z  IS  ROTATED  120  DEGREES 
COUNTER-CLOCKWISE  IF  ITS  IMAGINARY  PART  IS  NEGATIVE;  OTHER¬ 
WISE  Z  IS  ROTATED  120  DEGREES  CLOCKWISE.  THE  WRONSKIAN  OF  THE 
NUMERICALLY  INDEPENDENT  AIRY  FUNCTIONS  IS  STABLE  THROUGHOUT  THE 
COMPLEX  PLANE. 


GIVEN  Z.  IF  INDEP  <  0,  CLAIRY  RETURNS  THE  FIRST  AND  2ND  KIND 
AIRY  FUNCTIONS: 

(  AI,  A I P  }  =  EXP  ( ZT  AM)  *  (  AHZ),  AHZ)  ) 

$ 

(  BI  ,  BIP  }  =  EXP ( -ZTAM)  •  (  BI(Z),  BI (Z)  ).  IF  REAL(ZTAM)  >  =  0 

=  EXP(ZTAM)  *  (  BHZ),  BUZ)  ),  IF  REAL(ZTAM)  <  0 

SINCE  AHZ)  -  -/  +  1BHZ)  AND  A1(Z)  -  -/+  1BHZ)  IF  60  <  ABS C ARG  Z 
<  180.  THE  WRONSKIAN  OF  THE  AIRY  FUNCTIONS  OF  THE  FIRST  AND 

SECOND  KIND  IS  NOT  STABLE  THROUGHOUT  THE  COMPLEX  PLANE. 

***************************************************************** 

SUBROUTINE  CLAIRY (Z ,  AI,  BI ,  AIP,  BIP.  ZTAM,  INDEP  ) 

IMPLICIT  REAL*8  (  A-H,  O-Z  ) 

COMPLEX* 16  ZS.  ZSS.  ZTA ,  SCI,  SC2,  SD1,  SD2,  SQM1,  TSQM1 , 

+  HEPI06,  CHEP06 ,  EPI06,  CEPI06,  E2PI03,  CE2P03,  TEPI06,  TCPI06 

+  Tl,  T2 ,  CK,  ZSQ.  F,  G,  Q.  R,  ZB,  ZERO.  ZO.  REPI06.  E2ZTAM, 

+  Z.  AI .  BI,  AIP,  BIP,  ZTAM 
DIMENSION  C  (  40  )  ,  DUO) 

DATA  RDEG/ . 1745329251994329D-1/ ,  PI 04/ . 785398 1633974483D0/ 

DATA  RC2RPI / . 28209479 1773878 14D0/ ,  RPI / 1 . 7724538509055 16D0/ 

DATA  SQR3/  1 . 7 3 205 080 7 5688 7 7 29D0/  ,  ZERO  /  (  0 .  DO  ,  0  .  DO)  / 

DATA  SQM1/CO.ODO,  l.ODO)/.  TSQM1/ (0 . ODO ,  2.0D0)/ 

DATA  SN60 , CS60/ . 8 660 2540 37 644 3864 7D0 ,  0.5D0/, 

+  EPI 06/ C . 6660 25403 7 844 38 64 7 DO ,  O.5D0)/. 

+  CEPI06/C. 8 6602 5403 7 8 44 3 8 64 7D0 ,  -0.5D0)/, 

+  TEPI06/C 1. 7 3 20508075688 7 729D0 ,  l.ODO)/, 

+  TCPI06/( 1. 7 320508075688 7 7 29D0 ,  -l.ODO)/, 

*  E2PI03/(-0.5D0,  . 8 66025403 784438 64 7D0)/ , 

+  CE2P03/ ( -0 . 5D0 ,  - . 866025403784438647D0) / . 

*  HEPI06/C . 4330 1270 1B922 193233D0 ,  0.25D0)/, 

*  CHEP06/C .4330 1270 18922 193233D0,  -0.25D0)/ 

DATA  C/  1 . DO ,  0 . 69444444444444444D-0 1 . 

+  0.5347 2222222222222D+ 00,  0 . 1023 148 148 148 148 1D+01 , 

+  0. 15173611111111111D+01,  0.-20138888888888888D+01 , 

+  0.25115740740740740D+01,  0 . 30099206349206349D+01 , 
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+  0.350868055 55555555D+01, 
+  0 . 45069444444444444D+0 1 , 
+  0. 55057870370370370 D+01, 
+  0. ■65049603 174603 174D+01. 
+  0. 750434027777777770+01. 
+  0. 850385802469 13580D+O1. 
+  0.95034722222222223D+01, 
+  0. 10503156565656565D+02, 
+  0. 11502893518518518D+02, 
+  0. 12502670940170940D+02, 
+  0. 13502480158730158D+02. 
+  0. 14502314814814814D+02. 
+  0. 155021701388888880+02. 
+  0. 16502042483660130D+02. 
+  0. 1750 19 29 01 2345679 D+ 02, 
+  0. 18501827485380117D+02, 

DATA  0/  1 .  DO  , 

+  0.45 138888888888888 D+ 00. 
+  0. 14756944444444444D+01, 
+  0 . 24837962962962962D+0 1 . 
+  0 . 348784 72222222222D+01 . 
+  0 . 44902777777777777D+01 , 
+  0. 549189814814814820+01, 
+  0.64930555555555555D+01, 
+  0.7 49 39236111111111 D+01, 
+  0.8 4945987 65 4320987 D+01, 
+  0.949513888888888880+01. 
+  0 . 104955808080808080+02 . 
+  0.1 1495949074074074D+02 . 
+  0. 12496260683760683D+02, 
+  0. 13496527777777777D+02, 
+  0 . 14496759259259259D+02 , 
+  0. 15496961805555555D+02, 
+  0. 16497 1405228758 17D+02 , 
+  0. 174972993827 16049D+02 , 
+■  0. 18497441520467836D+02. 


0. 40077 1604938 27 160D+01. 
0.50063131313131313D+01, 
0.600534188034188030+01, 
0 . 700462962962962960+0 1 , 
0.800408496732026140+01. 
0 . 90036549707602339D+0 1 . 
0. 10003306878306878D+02, 
0. 1 10030 1932367 1497 D+ 02, 
0. 12002777777777777D+02, 
0. 13002572016460905D+02. 
0. 14002 39 4636015325D+02 , 
0. 15002240 143369175D+02 , 
0 . 16002 104377 104377D+02 , 
0.170019841269841260+02, 
0. 18001876876876877D+02. 
0.190017806267806260+02/ 

-0.97222222222222222D-01, 

0 . 96759259259259259D+00 . 
0. 198055555555555550+01 , 
0.298611111 1111 1111D+01. 
0.3989 197530864 19750+01. 
0. 499 1 1616 1616 1616 1D+0 1 , 
0.599252 136752 13675D+01, 
0. 69935 185 185 185 185D+0 1 , 
0.799428 10457516340D+01. 
0 . 89948830409356725D+0 1 . 
0.999537037037037040+01, 
0. 109957729468599030+02. 
0. 119961111111111110+02, 
0 . 12996399 176954732D+02 , 
0 . 139966475095785440+02 . 
0. 14996863799283154D+02, 
0. 159970538720538720+02, 
0. 16997222222222222D+02. 
0. 17997372372372372 D+ 02, 
0. 18997507 122507 122D+02/ 


ZTAM  =  ZERO 

IF  (CDABS  C  Z)  .GT.  5.)  GO  TO  200 

POWER  SERIES  EXPANSION  OF  AIRY  FUNCTION  FOR  ABS(Z).LE.5 
T 1  =  DCMPLX ( .35502B053BB78 17239D0 , 0 . DO ) 

T2  =  DCMPLX ( . 2588 1940379280679DO , 0 . DO ) 

F  =  T1 

R  =  T2 

T2  =  T2*Z 

G  =  T2 

ZSO  =  Z*Z 

0  =  ZERD 

XN  =  1 . DO 


100  XN  =  XN+l.DO 

T 1  =  Tl+ZSO/XN 

Q  =  Q+Tl 

XN  =  XN+l.DO 

T 1  =  Tl+Z/XN 

CK  =  F 

F  =  F+Tl 
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T2  =  T2*ZSQ/XN 

R  =  R+T2 


XN  =  XN+l.DO 

T2  =  T2*Z/XN 

G  =  G+T2 

C 

IF  (  DREAL(CK) .NE.  DRE A L ( F ) . OR . D I MAG ( CK ) . NE . DI MAG ( F ) )  GO  TO  100 
C 

A I  =  F-G 
A I  P  =  Q-R 
BI  =  SQR3* (F+G) 

BIP  =  SQR3* (Q+R) 

I F (  INDEP  .GT.  0  )  THEN  !  SCALED  INDEPENDENT  BIRY,  AIRY 
I F (  DIMAG(Z)  .GE.  O.ODO  )  THEN 

BI  =  HEPI06  *  (  BI  -  SQM1  «  AI  ) 

BIP  =  HEPI06  *  (  BIP  -  SQM1  •  AIP  ) 

ELSE 

BI  =  CHEP06  *  (  BI  +  SQM1  *  AI  ) 

BIP  =  CHEP06  *  (  BIP  +  SQM1  *  AIP  ) 

END  IF 
END  IF 
RETURN 
L 

C  ASYMPTOTIC  EXPANSION  OF  AIRY  FUNCTION  FOR  LARGE  Z. 

C  GET  ABS (  ARG  Z  ) 

200  TH  =  DABS (  DATAN2(  DIMAG(Z) ,  DREAL(Z)  )  /  RDEG  ) 

T2  =  E2PI03 

IF(DIMAG(Z)  .GE.  O.ODO)  T2  =  CE2P03 

C 

IF(TH  .GT.  150. ODO)  THEN 

ZO  =  Z  *  T2  !  ROTATE  Z 

ZS  =  CDSQRT ( ZO) 

ZSS  =  CDSQRT ( ZS) 

ZTA  =  ZO  *  ZS  /  1.5D0 
ZTAM  =  -ZTA 
E2ZTAM  =  ZERO 

IF  (  DREAL (ZTA)  .LE.  43. DO  )  E2ZTAM  =  CDEXP(2.D0  *  ZTAM) 

CALL  SUMMER ( ZTA ,  SCI,  SC2,  C,  0) 

CALL  SUMMER (ZTA ,  SD1.  SD2 ,  D.  1) 

BI  =  RC2RPI  /  ZSS  *  SCI  !  EXP(ZTA ) A I (ZO ) 

BIP  =  -RC2RPI  *  ZSS  *  SD1  «  T2  !  EXP(ZTA)A1p(Z0)T2 
I F (DIMAG(Z )  .GE.  O.ODO)  THEN 

AI  =  CHEP06  *  (SC2/ (RPI *ZSS)  ♦  ( TSQM1*E2ZT AM) *BI ) 

AIP  =  CHEP06  *  (SD2/RPI*ZSS*T2  ♦  (TSQM1*E2ZTAM) *BI P) 

ELSE 

AI  =  HEPI06  *  (SC2/(RPI*ZSS)  -  (TSQM1*E2ZTAM) *BI ) 

AIP  =  HEPI06  *  (SD2/RPI*ZSS*T2  -  (TSQM1*E2ZTAM) *BI P) 

END  IF 

ELSE  !  FOR  ABS(  ARG  Z)  .LE.  150  DEGREES 

C  USE  NBS  10.4.59,  10.4.61  FOR  AI ,  AIP 

ZS  =  CDSQRT (Z) 

ZSS  =  CDSQRT ( ZS ) 

ZTA  =  Z*ZS/1.5D0 

CALL  SUMMER  (ZTA,  SCI,  SC2,  C,  0) 

ZTAM  =  ZTA 

AI  =  RC2RPI/ZSS*SC1 

CALL  SUMMER  (ZTA,  SD1,  SD2 ,  D,  1) 

AIP  =  -RC2RPI*ZSS*SD1 
REPI06  =  EPI06 

O 

CM 
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I F (  DIMAG(Z)  .GE.  O.ODO  )  REPI06  =  CEPI06 
BIP  =  -RC2RPI  *  SD2  «  ZSS  *  REPI06  «  T2 
BI  =  RC2RPI  *  SC2  /  (ZSS  *  REPI06) 

END  IF 
C 

I F (  INDEP  .LT.  0)  THEN  !  SCALED  AIRY  OF  2ND  KIND 
E2ZTAM  =  O.ODO 

I F (  DREAL(ZTAM)  .GE.  O.ODO  )  THEN  !  ABS(  ARG  Z  )  <=  60 

I F (  DREAL(ZTAM)  .LE.  43. DO  )  E2ZTAM  =  CDEXP(-2.D0  »  ZT  AM ) 
I F (  DIMAG(Z)  .GE.  O.ODO  )  THEN 

BI  =  SC 2  /  (ZSS*RPI )  +  ( SQM1  «  E2ZTAM)  «  AI 
BIP  =  SD2  •  ZSS  /  RPI  +  ( SQM  1  «  E2ZTAM)  •  AIP 
ELSE 

BI  =  SC2  /  (ZSS*RPI )  -  ( SQM1  *  E2ZTAM)  *  AI 
BIP  =  SD2  *  ZSS  /  RPI  -  ( SQM  1  «  E2ZTAM)  «  AIP 
END  IF 

ELSE  !  60  <  ABS (  ARG  Z  )  <=  180 

I F (  DREAL(ZTAM)  .GE.  -43. DO)  E2ZTAM  =  CDEXP(2.D0  *  ZTAM) 

I F (  DIMAG(Z)  .GE.  O.ODO  )  THEN 

BI  =  (TCPI06  «  E2ZTAM)  *  BI  +  SQM1  •  AI 
BIP  =  (TCPI06  *  E2ZTAM)  *BIP  +  SQM1  *  AIP 
ELSE 

BI  =  (TEPI06  «  E2ZTAM)  «  BI  -  SQM1  «  AI 
BIP  =  (TEPI06  «  E2ZTAM)  «  BIP  -  SQM1  «  AIP 
END  IF 
END  IF 
END  IF 
RETURN 
END 
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SUBROUTINE  SUMMER  ( ZTA ,  SI.  S2,  CF ,  NGATE ) 

DOUBLE  PRECISION  CF 
REAL  MSTAR 

COMPLEX* 16  ZTA.VZTA.G.S1.S2 
DIMENSION  CF ( AO ) 

DATA  MSTAR/19.35/ 

MSTAR  IS  THE  LIMIT  OF  THE  SEQUENCE: 

M ( J )  =  -LN(  M(J-l)  )  /  4  +  M( 0 ) 

C  M( 0 )  =  -LN(  7.2  U  /  SQR(PI)  )  /  2  ) 
WHERE  U  IS  THE  COMPUTER'S  UNIT  ROUND¬ 
OFF  ERROR  (U  =  2** (-60)  HERE). 

I F (NGATE . EQ . 0 ) THEN 

ABSZTA  =  COABS ( ZTA ) 

I F ( ABSZTA  .GT.  MSTARJTHEN 

FETCH  NUMBER  OF  TERMS  FROM  TABLE 
NTERMS  =  NTERMC  ABSZTA  ) 

END  IF 
END  IF 

I F ( ABSZTA  .LE.  MSTARJTHEN 
NTERMS  =  2 . 0*ABSZTA  +  2.0 

IFCSNGLCCF (NTERMS ) )  .GT.  ABSZTA)  NTERMS  =  NTERMS  -  1 
END  IF 

G  =  DCMPLXCCFC 1) ,0.000) 

51  =  G 

52  =  G 

VZT  A  =  G/ZTA 

DO  100  I  =  2,  NTERMS,  2 

G  =  G* (VZTA *CF ( I )  ) 

51  *  SI  -  G 

52  =  S2  +  G 

G  =  G*(VZTA*CF(I+1) ) 

51  =  SI  +  G 

52  =  S2  +  G 
100  CONTINUE 

IFC  MODCNTERMS.  2)  .NE.  0  )  RETURN 

51  =  SI  -  G 

52  =  S2  -  G 
RETURN 

END 
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FUNCTION  NTERM(ABSZTA) 

C  ***  USING  THE  KEY  ABSZTA ,  FUNCTION  NTERM  DOES  A  BINARY  SEARCH 
C  **«  OF  AZT A  TO  DETERMINE  THE  NUMBER  OF  TERMS  OF  THE  SERIES  TO 
C  *•*  SUM.  THE  TABLE  OF  AZTA  VALUES  WAS  CONSTRUCTED  FOR  U  = 

C  ***  2**(-60).  IF  THE  VALUE  OF  U  CHANGES  THEN  A  NEW  TABLE  OF 
C  *«*  AZTA  VALUES  SHOULD  BE  CONSTRUCTED  AS  DESCRIBED  IN  NUSC 
C  ***  TM  861032 . 

C  *** *********************************************************** 
INTEGER  NTERM 

DIMENSION  AZTA ( 3B ) ,  NTABLE (38 ) 

DATA  (AZTA(K).  NT  ABLE ( K) ,  K  =  1.  19)/ 

+  19.35  .39. 

+  19.36194.  38, 

+  19.38579,  37, 

+  19.42577,  36, 

+  19.48381,  35, 

+  19.5621  ,  34, 

+  19.66328.  33, 

+  19.79037,  32, 

+  19.94698,  31, 

+  20.1374  ,  30, 

+  20.36676,  29, 

+  20.64125,  28. 

+  20.96842,  27. 

+  21.35754,  26. 

+  21.8201  ,  25, 

+  22.37055,  24, 

+  23.02719,  23. 

+  23.81359,  22, 

+  24.76054,  21/ 

DATA  (AZTA(K),  NTABLE ( K) ,  K  =  20.  38)/ 

+  25.90883,  20. 

+  27.31371,  19, 

+  29.05135,  18, 

+  31.22937,  17, 

+  34.00401,  16, 

+  37.60926,  15. 

+  42.4084,  14, 

+  48.99014,  13, 

+  58.35861,  12, 

+  72.336  ,  11, 

+  94.49151,  10, 

+  132.5262  ,  9, 

+  205.2731  ,  8, 

+  367.2325  ,  7, 

+  818.2973,  6, 

+  2605.586  ,  5. 

+  15653.95  ,  4, 

+  342085.9  ,  3, 

+  342085.9  ,  3/ 

I  =  1 

'  J  =  37 

C  DO  WHILE  (  I  .LE.  J  ) 

20  I F  (  .NOT.  (I  .LE.  J)  )  GO  TO  70 

L  =  (  I  +  J  )  /  2 
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I F (  ABSZTA  .LT.  AZTA(L)  )  THEN 
J  =  L  -  1 
ELSE 

I  =  L  +  1 
END  IF 
GO  TO  20 
70  CONTINUE 

NTERM  =  NTABLE(  I  ) 

RETURN 

END 
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